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ABSTRACT 

We present a simulation of the long-term evolution of a Population III supernova remnant 
in a cosmological minihalo. Employing passive Lagrangian tracer particles, we investigate 
how chemical stratification and anisotropy in the explosion can affect the abundances of the 
first low-mass, metal-enriched stars. We find that reverse shock heating can leave the inner 
mass shells at entropies too high to cool, leading to carbon-enhancement in the re-collapsing 
gas. This hydrodynamic selection effect could explain the observed incidence of carbon- 
enhanced metal-poor (CEMP) stars at low metallicity. We further explore how anisotropic 
ejecta distributions, recently seen in direct numerical simulations of core-collapse explosions, 
may translate to abundances in metal-poor stars. We find that some of the observed scatter 
in the Population II abundance ratios can be explained by an incomplete mixing of super¬ 
nova ejecta, even in the case of only one contributing enrichment event. We demonstrate that 
the customary hypothesis of fully-mixed ejecta clearly fails if post-explosion hydrodynamics 
prefers the recycling of some nucleosynthetic products over others. Eurthermore, to fully ex¬ 
ploit the stellar-archaeological program of constraining the Pop III initial mass function from 
the observed Pop II abundances, considering these hydrodynamical transport effects is crucial. 
We discuss applications to the rich chemical structure of ultra-faint dwarf satellite galaxies, to 
be probed in unprecedented detail with upcoming spectroscopic surveys. 

Key words: dark ages, reionization, first stars — galaxies: dwarf — galaxies: formation — 
methods: numerical — stars: abundances — stars: Population II 


1 INTRODUCTION 

The first generation of stars, the so-called Population III (Pop III) 
stars, formed from pure hydrogen, helium, and lithium (for re¬ 
view, see Bromm 2013, and references therein). Nucleosynthesis 
in these stars and their explosions delivered the first metals to en¬ 
rich the Universe. Then, second-generation Population II (Pop II) 
stars formed from sufficiently metal-enriched gas elements. The 
first generation of low-mass, long-lived Pop II stars, the relics of 
which can be found in the Galactic stellar halo and ultra-faint dwarf 
satellite galaxies (UFDs), encodes a crucial empirical record of 
Pop III nucleosynthesis (Frebel & Norris 2015). The principal vari¬ 
able controlling the stellar nucleosynthetic output is the progenitor 
mass. Simulations published over the past few years indicate that 
two processes can limit the masses of Pop III stars: radiative sup¬ 
pression of protostellar accretion (McKee & Tan 2008; Hosokawa 
et al. 2011; Stacy et al. 2012; Hirano et al. 2014; Susa et al. 2014) 
and gravitational fragmentation in protostellar disks (Stacy et al. 
2010; Clark et al. 2011a,b; Greif et al. 2011, 2012). An upper limit 
on the production rate of very-low-mass, long-lived Pop III stars 


can be placed empirically, by surveying for metal-free stars in the 
Galaxy (Tumlinson 2006; Hartwig et al. 2015). The failure to find 
Pop III stars in the nearby Universe suggests that the characteristic 
stellar mass was > O.SM©. Attempts have been made to specu¬ 
late about the Pop III initial mass function (IMF). On the basis of 
cosmological simulations and still uncertain subgrid prescriptions, 
Susa et al. (2014) estimated that the stars had masses from 1 to 
300 M©, that the IMF peaked at 45M©, and that there were on 
average 2.5 Pop III stars per cosmic minihalo. 

There are several ways in which the intuition gained in the 
study of metal-enriched star formation and evolution in the nearby 
Universe fails in metal-free stars. Pop III stars might have been 
more rapid rotators, as angular momentum loss through line-driven 
winds was not efficient in metal-free stellar atmospheres (Stacy et 
al. 2011).^ However, rapid rotation drives internal mixing of hy- 

^ It is premature to speculate about differences in magnetically-mediated 
angular momentum loss because the degree of magnetization of Pop III stars 
remains uncertain. 
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drostatically synthesized elements and this can dredge light ele¬ 
ments up into the atmosphere (Ekstrom et al. 2008; Yoon et al. 
2012). Also, stars lacking metals in their atmospheres have high 
enough surface temperatures ~ 10^ K to be strong ionizing pho¬ 
ton emitters (Schaerer 2002). The associated HII regions are ini¬ 
tially highly overpressured near the star. The ionized gas pressure 
reduces the density of gas into which the Pop III star would explode 
to ~ 0.1 — 1 cm“^. However it has recently become apparent that 
dense circumstellar gas clumps may avoid photoionization and re¬ 
main cold and dense until the supernova blastwaves arrive (Abel 
et al. 2007; Bland-Hawthorn et al. 2011; Ritter et al. 2012, 2015; 
Webster et al. 2014; Smith et al. 2015). These clouds modify the 
mechanics of blastwave-circumstellar medium interaction, for ex¬ 
ample by securing a reservoir of metal-free gas close to the site in 
which Pop II stars would eventually form. 

The Pop III explosion type and nucleosynthetic yields are sen¬ 
sitive to the mass of the compact remnant produced in the core 
collapse: the nucleosynthetic products collapsing into or accreted 
onto the remnant are removed from the yields. Gross theoretical 
uncertainties remain, especially with respect to the mass of the 
black hole produced by fallback and the metal fraction lost in the 
black hole. In non-rotating stars, the general expectation is of a 
core collapse supernova with a neutron star remnant for progeni¬ 
tor masses lOM© < M < 25Mq and a black hole remnant for 
25Mo < M < IdOM©. At still higher masses, a pair instability 
supernova (PISN) leaving no compact remnant or a direct collapse 
into a black hole are the likely outcomes (Nomoto et al. 2013). 
Stellar rotation tends to reduce the minimum mass for PISN events 
(Chatzopoulos & Wheeler 2012; Yoon et al. 2012). Estimates of 
the nucleosynthetic yields of Pop III supernovae as a function of 
progenitor mass are sensitive to additional theoretically uncertain 
parameters: the degree of rotation, mixing, and mass loss, the me¬ 
chanics of fallback (in particular as approximated with a simple 
mass cut), and explosion energy (Woosley & Weaver 1995; Heger 
& Woosley 2010; Takahashi et al. 2014). 

Particularly problematic are any conclusions derived from 
one-dimensional simulations of core-collapse supemovae. Multi¬ 
dimensional processes such as convection, instability, and mixing 
can be included in one-dimensional simulations with sub-grid pre¬ 
scriptions, but the dependence of the amplitude of these effects on 
progenitor parameters and any potential intrinsic stochasticity are 
poorly constrained. Eor example, explosively-synthesized nickel, 
the radioactive precursor to iron, can escape gravitational capture 
in the compact object by getting squirted ahead of the lighter ele¬ 
ments in collimated, Rayleigh-Taylor-instability-amplified jets (or 
“fingers”). This, of course, has pivotal implications for chemical 
enrichment. Under which specific conditions are fingers expected 
is the subject of recent, rapidly developing research (Hammer et al. 
2010; Wongwathanarat et al. 2013, 2015; Utrobin et al. 2014). 

At a similar stage of maturity as the study of metal-free stars, 
the narrative of the formation of the first low-mass, metal-enriched 
stars is still being written. The enhanced cooling that is required to 
facilitate gravitational fragmentation down to < 1 M© is provided 
by dust grains and dust-gas thermal coupling (Omukai et al. 2005; 
Schneider et al. 2006; Safranek-Shrader et al. 2014b, 2015). Dust is 
produced in the Pop III supemovae (Schneider et al. 2006; Nozawa 
et al. 2007; Cherchneff & Dwek 2010; Chiaki et al. 2015) or in red 
supergiant stellar winds (Nozawa et al. 2014). More speculatively, 
dust grows in gas phase in the collapse of metal-enriched clouds 
that would form Pop II stars (Chiaki et al. 2014). 

It has recently become possible to simulate the formation of 
individual metal enriched stars from first principles, from the ini¬ 


tial conditions imprinted in the Big Bang and by tracking all the 
relevant dynamical, chemical, and radiative processes (Safranek- 
Shrader et al. 2014b, 2015). These challenging simulations have not 
yet converged to making a definite prediction for, say, the shape of 
the Pop II stellar IMP, even in very-high-redshift systems in which 
the separation of the halo virial and star-forming time scales is the 
smallest and least computationally prohibitive. We are still some 
ways from a predictive theory of star formation across the halo 
mass scales and redshifts at which first-generation metal-enriched 
star formation is expected to take place. 

All these theoretical bottlenecks motivate the aspiration 
that the nucleosynthetic yields of Pop III stars can be reverse- 
engineered from chemical abundances measured (e.g., Cayrel et al. 
2004) in low-mass, long-lived stars in what can be called “stellar 
archaeology”. Because the most metal-poor stars may have formed 
from gas enriched by one or a small number of Pop III supemovae, 
measured stellar chemical abundances can be used to test theo¬ 
ries of Pop III nucleosynthesis (Prebel et al. 2005; Prebel & Norris 
2015). Theoretical predictions of isotopic yields are sensitive to pa¬ 
rameters not currently determined from first principles. Therefore 
the alternative approach of using the observed abundance ratios in 
the most metal-poor stars to identify classes of contributing super¬ 
novae in the chemical abundance space may be more promising 
(Ting et al. 2012; Tominaga et al. 2014). Eor instance, Ishigaki et al. 
(2014) have shown that the most iron-poor stars can be explained 
as having been enriched by low-energy Pop III core-collapse ex¬ 
plosions in which most of the iron was lost in the compact rem¬ 
nant. Indeed, faint supemovae with significant fallback seem to be 
necessary to explain the abundances of the most metal-poor stars 
(Iwamoto et al. 2005). The hypothetical ultra-energetic, magnetic- 
outfiow-powered “hypernovae” may explain abundance ratios in 
extremely metal poor stars (Izutani & Umeda 2010). The effects 
of rotation in Pop III stars may affect the C, N, and O abundances, 
but despite the expectation that Pop III stars were rapid rotators, 
very few nucleosynthesis calculations take rotation into account 
(Ekstrom et al. 2008; Joggerst et al. 2009, 2010). 

Metal-poor stars that could have been enriched by Pop III stars 
exist in the inner and outer stellar halo of the Milky Way, in old, 
metal-poor globular clusters, and in dwarf satellite galaxies. Espe¬ 
cially promising for stellar archaeology are the UEDs (Kirby et al. 
2008; Norris et al. 2010; Simon et al. 2010; Brown et al. 2012; Var¬ 
gas et al. 2013; Erebel et al. 2014). The inner and outer Galactic 
halo have peak metallicities of [Ee/H] = —1.6 and —2.2, respec¬ 
tively (York et al. 2000; Gunn et al. 2006; Carollo et al. 2008).^ 
The fraction of stars that are carbon enhanced and yet metal poor 
(CEMP, defined as [C/Ee] > +1.0 and [Ee/H] < —1.0) increases 
as the metallicity decreases; thus the outer halo has a higher CEMP 
fraction (Yong et al. 2013). Stars in UEDs have metallicity and 
abundance ratio distributions similar to those of the halo stars but 
with a few key differences (Corlies et al. 2013). The abundances 
of neutron capture elements in UEDs are lower than in halo stars 
(Erebel et al. 2010) while a-element abundances are more en¬ 
hanced in UEDs (Kirby et al. 2011). 

In addition to metal poor stars in the local group, the earliest 
stages of chemical evolution can be studied in intergalactic absorp¬ 
tion line systems such as the damped Lyman-a (DLA) systems on 
quasar (Cooke et al. 201 la,b, 2015; Becker et al. 2012) and gamma- 


^ We employ the standard notation: [A/5] = log(A/B) —log(A©/B©) 
where A and B are mass-fractions of isotopes A and B while A© and Bq 
are the corresponding solar mass fractions. 
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ray burst (Lamb & Reichart 2000; Wang et al. 2012; Cucchiara et 
al. 2015; Ma et al. 2015) sightlines. There has been an attempt to 
constrain the Pop III IMF by measuring the relative abundances of 
various elements in a sample of ~ 10 DLAs (Kulkarni et al. 2013). 
By comparing the abundance ratios in DLAs with those of halo and 
dwarf galaxy stars, we can learn how the systems that polluted the 
DLAs are related to the relic stellar systems in the nearby Universe 
(Salvadori & Ferrara 2012; Webster et al. 2015). 

Various explanations have been offered for the observed abun¬ 
dance patterns in metal poor stars and, in particular, the metallicity 
distribution function and the scatter in abundance ratios as a func¬ 
tion of [Fe/H]. Stars with the lowest metallicities may have formed 
from gas enriched by only one or possibly a few Pop III supernovae 
(Audouze & Silk 1995; Feltzing et al. 2009). The halo-to-halo vari¬ 
ation of Pop III properties (Hirano et al. 2014) could explain the rel¬ 
atively large scatter in abundance ratios that has been observed in 
the lowest metallicity stars (Frebel & Bromm 2012). The detailed 
abundances can be explained if most of the Pop III stars produced 
core collapse supernovae and only a few exploded as PISNe (Aoki 
et al. 2014). 

The decrease of the carbon-enhanced fraction with metallic¬ 
ity can be explained if more massive Pop III stars produce PISNe 
that synthesize more iron but also remove more gas from the halo 
(through photoionization and supernova blasting), thus making it 
less likely that iron-rich Pop II stars will subsequently form in that 
halo (Cooke & Madau 2014). The existence of carbon-normal and 
carbon-enhanced metal-poor stars has also been interpreted as evi¬ 
dence for two distinct cooling channels in early metal-enriched star 
forming systems. In this scenario, CEMP stars form from gas that 
is rich in C and O and first cools via atomic fine structure transi¬ 
tions, while carbon-normal stars form via non-carbon (e.g., silicate) 
dust-mediated cooling (Norris et al. 2013; Ji et al. 2014). The de¬ 
crease in abundance scatter with metallicity is then a consequence 
of an increase of the number of contributing nucleosynthetic events 
(Argast et al. 2000). However it remains unclear if the hypothetical 
dust-free, C- and 0-enriched gas can indeed efficiently fragment 
into low-mass protostellar cores. 

This is our third paper in a sequence (Ritter et al. 2012, 
2015) investigating the small-scale post-supernova evolution af¬ 
fecting interpretations of this record. Recognizing that the charac¬ 
ter of the Pop Ill-to-II transition can be sensitive to the poorly- 
explored small-scale hydrodynamic effects in supernova remnant 
(SNR) evolution, Ritter et al. (2012) followed the metal-enriched 
ejecta from a single core-collapse-type Pop III explosion in a 
10 ® M 0 minihalo—starting from the early, kinetic-energy domi¬ 
nated free expansion phase—for about 40 Myr after the explosion. 
They found that the metals were inhomogeneously dispersed in 
the host halo and its immediate neighborhood. The metal-enriched 
gas returning into the halo center became diluted with primordial 
gas streaming from the cosmic web only in a dense, numerically- 
unresolved gaseous core at the center of the simulation. 

Proceeding to a higher spatial resolution, Ritter et al. (2015) 
showed that turbulence excited by gravitational infall in a col¬ 
lapsing gaseous core is sufficient to homogenize chemical abun¬ 
dances in the highest-density gas. Feng & Krumholz (2014) showed 
that turbulent homogenization happens in collapsing, star-cluster- 
forming clouds in the present-day Galactic disk (see, also, Bland- 
Hawthorn et al. 2010). Interestingly, however, Ritter et al. (2015) 
also found that the high-density homogenization did not uniformly 
sample the abundance profiles initially injected by supernovae. 
There seemed to be a trend in which the yields that were raised 
to higher entropies in blastwave shocks were deficient in the fall¬ 


back and recycling into subsequent stellar generations. Ritter et al. 
suggested that this trend arose from a positive correlation between 
the cooling time in the supemova-ejecta-enriched gas and the spe¬ 
cific entropy to which the reverse supernova shock has raised the 
ejecta that would enrich the gas. 

Here, we continue analysis of one of the two simulations of 
Ritter et al. (2015) to find that abundance ratios in the ejecta of a 
Pop III supernova need not be preserved in star-forming clouds. A 
non-trivial, stochastic relation between yields and enrichment has 
to be taken into account when comparing the abundances measured 
in metal-poor stars with the predicted yields of supernovae. The 
mismatch arises from a combination of two hydrodynamical ef¬ 
fects, one related to the realization that the yields can be anisotrop- 
ically distributed in the ejecta, and the other that the reverse shock 
raises the inner ejecta to entropies too high to cool and be recycled 
into new stars. Whereas Ritter et al. (2015) focused on the case 
of multiple supernovae, we here treat the case of inhomogeneous 
and anisotropic mixing biases from a single Pop III explosion. We 
thus derive, as it were, the irreducible degree of such biases, and 
investigate how they translate into a minimal degree of scatter in 
the empirical abundance record. 

The rest of the paper is organized as follows. In Section 2 we 
describe our numerical methodology, including in particular detail 
left out in Ritter et al. (2015). In Section 3 we discuss morphologi¬ 
cal evolution of the SNR, ejecta dispersal trends, chemical variation 
in diffuse and dense gas, and explore implications for identifying 
chemical abundance biases and anomalies in observations of metal- 
poor stars. In Section 4 we briefly refiect on the broader state of 
stellar archaeology. 

2 METHODOLOGY 

2.1 Initial conditions, dark matter, and hydrodynamics 

The initial conditions were mapped from a snapshot of the 
gravitational-hydrodynamic simulation reported in Ritter et al. 
(2012). Briefiy, the parent simulation had been carried out in a 
periodic cosmological box of comoving volume IMpc®. Gaus¬ 
sian perturbations in the parent simulation were initialized at red- 
shift of z — 146 with the code GRAPHIC2 of Bertschinger 
(2001). Chemical composition of the gas was evolved with a non¬ 
equilibrium chemical network for the primordial chemical species 
as in Safranek-Shrader et al. (2012). At redshift z = 19.7, a 
10 ® M© virialized halo formed and near the center of the halo the 
gas density reached n ~ 10® cm“®. At this time, we inserted a star 
particle at the density maximum. The star particle accelerated un¬ 
der the combined gravitational force of the gas and dark matter but 
was otherwise collisionless. 

To initialize the present simulation, we took the cosmological 
simulation at the instance of star particle insertion and excised a 
cubical region of physical width 1 kpc centered on the particle. Hy¬ 
drodynamic variables in the region were mapped onto an isolated, 
non-expanding cubical box of the same size as the excised region. 
The center-of-mass velocity in the mapped domain was subtracted 
from the gas velocity in each cell. The density, temperature, and 
pressure were mapped directly from the parent grid. We did not 
map the dark matter particles. Instead, we modeled the dark matter 
mass distribution with a time-dependent power-law density profile 
that was spherically symmetric around the box center. From the in¬ 
sertion of the star particle (the beginning of the simulation) until 
the insertion of the SNR, the dark matter density profile was held 
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constant 

PDM(r) = 940rpc^‘^^^Mopc“^ (1) 

where r^c is the distanee from the center of the halo in parsecs. 
This profile was obtained as a direet spherically-symmetric power- 
law fit to the dark matter density in the parent simulation at the 
point of star particle insertion. To model dark matter halo growth 
from 1 Myr after SNR insertion onward, we implemented a time 
dependence through a scaling calibrated against the halo growth 
rate in the Ritter et al. (2015) continuation of the simulation over 
200 Myr after supernova insertion 

PDM(r, t) = f{t) g{t)~^pYiM(,g{t)~^r), (2) 

where f{t) = 1.029*“^', p(t) = 1.012‘“y% and ^Myr is time after 
supernova insertion in Myr. 

The simulation was carried out using the adaptive mesh re¬ 
finement (AMR) hydrodynamics code FLASH, version 4.0 (Fryx- 
ell et al. 2000; Dubey et al. 2008). We used the directionally-split 
piecewise-parabolic method in three dimensions to solve the hy¬ 
drodynamics equations. The gravitational potential was computed 
by adding the analytical dark matter density profile in Equations 
(1) and (2) to the gas density on the AMR grid and then solving 
the Poisson equation with the FLASH multigrid iteration scheme 
(Ricker 2008) with isolated boundary conditions. The simulation 
was not cosmological in the sense that it did not include universal 
expansion and active dark matter particles. The size of the mapped 
box was just large enough to avoid corruption of the solution at the 
box center by boundary artifacts over the 60 Myr of evolution. 


2.2 Pop III stellar model, photoionization, and cooling 

We used a single collisionless star particle to represent the position 
of the Pop III star cluster that we assumed had formed at the center 
of the halo. The star particle acted as an isotropie source of ionizing 
radiation with luminosity Q(H) = 6 x 10^^ ionizing photons s“^ 
over a stellar lifetime of 3 Myr. The particle was placed at the gas 
density maximum at the beginning of the simulation. For refer¬ 
ence, if Q(H) was to be emitted by single metal-free star, its mass 
would have been ~ 70 M© (Schaerer 2002). We consistently set 
our ejeeta mass to be lower than this, 40 M©, to allow for mass 
loss to a compact remnant and for additional, longer-living Pop III 
stars that could have formed in the primary star’s protostellar disk 
and eontributed their own ionizing radiation (see Seetion 2.3). 

We assumed photoionization equilibrium in cells exposed to 
ionizing radiation from the star and collisional ionization equilib¬ 
rium in all other eells. For photoionization equilibrium we used the 
code CLOUDY (Ferland et al. 2013) to tabulate the temperature and 
ionization state as a function of the ionization parameter, ^ ^ F/n, 
where F is the ionizing photon flux and n is the gas number density. 
For collisional ionization equilibrium, we used CLOUDY to tabulate 
the ionization state and cooling rate as a function of density, tem¬ 
perature, and metallicity. The tabulated data was interpolated at run 
time. Cooling was operator-split from the hydrodynamic update. 
Where necessary, cooling was sub-cycled on time steps equal to 
one-tenth of the local cooling time. 

As a simplification, we did not include molecules and dust 
grains in the equilibrium calculations.The tracking of molecule for¬ 
mation would have required the integration of a non-equilibrium 
ehemical rate equation, something we did in the parent cosmo- 
logieal simulation but omitted here for the sake of eomputational 
efficiency. The absence of H 2 molecules in the recollapsing rem¬ 
nant might have artificially delayed gas collapse past densities 


~ 10cm“^. Had H 2 or dust been present, the runaway gravita¬ 
tional collapse could have potentially set in much earlier. 

The details of ionizing radiative transfer were the same as in 
Ritter et al. (2015). Briefiy, the eartesian AMR grid was mapped 
onto a spherical grid defined with the HEALPix scheme (Gorski 
et al. 2005). Along eaeh of the ~ 3000 angular pixels and out to 
eaeh pixel-specifie Stromgren radius, we transported the ionizing 
flux between bins and then mapped it back onto the AMR grid. In 
the CLOUDY calculations, the souree was assumed to have a Planck 
spectral energy distribution at Tgff = 10^ K. The method is not pre¬ 
cisely photon-conserving but was sufficient to reproduce the correct 
hydrodynamic response of the gas to photoionization heating (ra¬ 
diation pressure was not computed). Figure 1 shows the effect of 
ionization on the density and temperature of the gas in the halo. 
When the star particle is inserted, the gas density at the eenter of 
the halo is ~ 60 em“^ and the temperature is ~ 1000 K. The tem¬ 
perature quickly rises to a central maximum of ^ 70,000 K and the 
pressure excess drives the density to gradually drop to ~ 0.2 cm“^ 
within the central few par sees after 3 Myr. 


2.3 Supernova remnant and ejecta tracking 


The ionizing source was turned off 3 Myr after star particle inser¬ 
tion and a SNR was inserted at the location of the star particle. The 
initial radius of the remnant of Rqn = 1.87 pc was chosen such 
that the remnant was in free expansion and the blastwave has not 
had the chance to sweep a circumstellar gas mass comparable to the 
ejecta mass. The SNR was inserted by modifying the gas density, 
temperature, absolute metallicity, and gas velocity in a sphere of 
radius i^sN = 1.87 pc. We followed the prescription of Whalen et 
al. (2008) to set the density and velocity to 


P = 
and 

V = 


Pcore 

Pcore 



L Rcore^ 
L ^ Rcore') 


^max / 


rsN), 


(3) 

(4) 


where rsN denotes the center of the SNR, Wmax = 10® cm s ^, and 
(3 — The SNR core density pcore was computed with 


Pcore 


3Mej / I3\ 

V 3 / 


X 


^ (R core 

3 \ Rsn 


1.42 X 10 gem 




(5) 


and the core radius with 


Rsn I f0{l3 — 5)F^kin 
^max y 3(^ - 3)Mej 

9.62x10^^ cm, (6) 

= 10^^ erg is the kinetic energy of the explosion and 
Mej = 40 M© is the ejecta mass. The temperature was set to an 
arbitrary uniform initial value of 10^ K. Since we do not simulate 
molecule and dust formation, the temperature of the ejecta preced¬ 
ing reverse shock passage, as long as it is low enough (^ 10^ K) to 
ensure that kinetic energy dominates, is irrelevant for the physical 
outcome of the simulation. 

The instanee of SNR insertion eorresponds to a physical time 
of RsN/vma^ = 182 yr after core-collapse explosion. The total 
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Figure 1. Radial profiles of density and temperature at four times before and two times after SNR insertion. The profiles were spherically averaged around the 
star particle before insertion and around the insertion site after insertion. The time since the beginning of the simulation is indicated in the legend. The blue 
line labeled 3 Myr shows the time immediately before SNR insertion and the purple line also labeled 3 Myr shows the time immediately after SNR insertion. 
The brown line labeled 3.034 Myr shows the time when the reverse shock has first reached the center of the halo. Also shown as dotted lines are the initial 
core density pcore, the core radius Rcore, and the initial SNR radius as well as the virial radius i^vir and temperature Tvir of the halo. 


pre-explosion gas mass displaced by SNR insertion was 0.08 Mq, 
negligible compared to the ejecta mass. The small displaced mass 
justifies our neglect of the blast wave-circumstellar matter interac¬ 
tion in the period preceding insertion. 

We track the metal-bearing ejecta redundantly in two ways: 
with an Eulerian passive scalar variable representing the absolute 
metallicity, and with a large number of Lagrangian tracer parti¬ 
cles with initial distribution matching the ejecta density. Redundant 
tracking is required to control the artificial numerical diffusion to 
which Eulerian passive scalar advection is susceptible. This dif¬ 
fusion introduces spurious quasi-exponential low metallicity tails 
where metal-free gas is in contact with metal rich gas (see Ritter 
et al. 2012, 2015, for further discussion of numerical diffusion in 
passive scalar transport).^ 

We assumed that 4 M© of the ejecta was in metals. We com¬ 
pute the cooling rate and ionization state of the gas from the local 
passive scalar metallicity which for simplicity we treated as fully 
mixed within the initial remnant so that metallicity within the rem¬ 
nant was Z = 0.1 bZ© and abundances were a-enhanced 
([a/Fe] = 0.5) but otherwise solar (in contrast with the metal 
transport analysis in Section 3 where mixing and a-enhanced solar 
abundances were not assumed). 

At 1 kyr after the SNR was inserted in the simulation, 10^ 
passive tracer particles were also inserted to track the motion of 
the ejecta. The particles were advected with the gas velocity cubi- 
cally interpolated from the computational grid. Each particle repre¬ 
sented a mass of Mpart = 4 x 10“® M© of ejecta gas. The num¬ 
ber of particles placed in each AMR cell of volume AV was set 
to a random integer having the same mean as, and within unity of 


^ Obergaulinger et al. (2014) have recently employed the tracer particle 
technique to study ejecta hydrodynamics in the Vela Jr. SNR. 


Z p Ay/Mpart, where Z and p are the passive scalar metallicity 
and density in the cell. 

Each Lagrangian tracer particle was given a unique tag so that 
the particle’s location can be traced for the duration of the simu¬ 
lation. A spherical coordinate system (r, 0) was centered on the 

explosion with 0 being the inclination with respect to the positive 
z-axis and 0 the azimuth relative to the positive x-axis. The par¬ 
ticles were then classified by these coordinates. At the end of the 
simulation particles that were located in cells with gas densities ex¬ 
ceeding 10“^^ g cm“^ were treated as having the chance to having 
been incorporated in a Pop II star cluster that would form at the 
center of the halo if the hydrodynamical evolution were tracked to 
higher densities and spatial resolutions. Particle tagging allowed us 
to trace every particle that made it into the Pop II cluster to its initial 
location in the supernova ejecta. 

2.4 Adaptive refinement 

While the ionizing source was present in the simulation, the compu¬ 
tational block (8x8x8 computational cells) containing the source 
was forced to be at refinement level ^max = 7 and the cell size 
was Axmin = 1 kpc/2^™'"^^^ 2 pc. The other blocks were at 

the lowest refinement level possible subject to the FLASH constraint 
that adjacent computational blocks can differ in refinement level by 
at most one (| A^| ^ 1). At the end of the star’s life (tsN = 3 Myr), 
the hydrodynamic time step was reset to a small value, the ionizing 
source was turned off, and the block containing the star particle was 
refined so that the particle now found itself in a block at a maximum 
refinement level of ^max = 13 with cell size Aa^min ~ 0.03 pc. The 
region at maximum refinement was made just large enough to let 
the entire SNR be inserted at the same, uniform refinement level. 
As before, the blocks not required to be at the maximum refinement 
level were derefined as much as AMR constraints allowed. 
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Figure 2. The columns show the evolution of the remnant at three different times after insertion, where the time since the beginning of the simulation is shown 
in the bottom left corner of each plot (the SNR was inserted at 3 Myr). The rows show: density slices (top row), density-weighted gas density projections 
(f pndzj f pdz; second row), passive scalar metallicity slices (third row), and density-weighted passive scalar metallicity projections (f pZdzj f pdz; 
bottom row). Note that the metallicity color scale is truncated to display metallicities exceeding 0.01 Zq with the same color even where the local metallicity 
substantially exceeds this. 
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t - tsN (Myr) 


Figure 3. The median radii of Lagrangian tracer particles originally at radial 
mass coordinates <0.01 Mej (red lines), 0.1, 0.2, ..., 0.9 Mej (orange to 
green line), and > 0.99 Mej (purple line) as a function of time from SNR 
insertion. Note the inversion in which inner shells overtake outer shells at 
late times. 

After the supernova remnant was inserted, the hydrodynamic 
time step was allowed to increase freely subject to the Courant- 
Friedrichs-Lewy criterion. The extent of the region forced to max¬ 
imum refinement was progressively enlarged such that the entire 
SNR as identified by the extent of the forward shock remained at 
maximum refinement. The maximum refinement level ^max was 
allowed to decrease subject to the requirement that the diameter of 
the SNR be resolved by at least 128 cells. This meant that the max¬ 
imum refinement changed to ^max = (12, 10, 9, 8 ) ait — tsN = 
(0.0044, 0.18, 1.2, 8.7) Myr, respectively, with ^max = 8 being 
the maximum refinement level at the end of the simulation. At late 
times t — tsN > 8.7 Myr when we focused on the gravitational 
fallback of metal-polluted gas to the center of the halo, we did not 
keep the entire remnant at the same refinement level, but instead 
refined the grid using the standard FLASH second-derivative-based 
refinement criterion applied to temperature, density, and metallic- 
ity. The Jeans length in the metal-enriched gas was resolved by at 
least 16 cells until the end of the simulation. 

3 RESULTS 

3.1 Morphological evolution of the supernova remnant 

Figure 2 shows slices and projections of density and passive scalar 
metallicity at three different times after SNR insertion. In the left 
column showing the state of the remnant at 0.9 Myr after the 
explosion, we see a thin snowplow shell enclosing a low-density 
(~ 10“^ — 10“^ cm“^) interior. The elliptical features exterior 
to the snowplow shell are the shock wave driven outward by HII 
region pressure. The top-left-to-bottom-right diagonal feature is a 
gaseous filament of the cosmic web. Metals remain confined to the 
SNR shell and its interior. 

At 35 Myr after explosion, hydrodynamic instability has al¬ 


ready rendered the SNR highly anisotropic. It has expanded to the 
box edges (> 500 pc) and has collapsed in the direction of the cos¬ 
mic web filament. The supernova ejecta and entrained gas have be¬ 
gun to stream back into the halo center. The densest gas is con¬ 
fined by the dark matter halo’s gravitational pull and is not self- 
gravitating. The peak gas density remains low, ~ 1 — 10 cm“^, in¬ 
sufficient to provide conditions for star formation. Metal-carrying 
ejecta are confined in a thin, crumpled sheet perpendicular to the 
cosmic web filament. The ejecta remain highly inhomogenous 
and undiluted: this may be less clear in the numerical-diffusion- 
susceptible passive scalar metallicity slice shown in Figure 2, but 
is visible in the Lagrangian tracer particles (not shown) that are 
distributed in thin sheets. 

At 57 Myr after explosion, the dark-matter-confined, turbu¬ 
lent, quasi-hydro static gas has reached density 10^ cm“^. The gas 
Jeans length has dropped to Aj ~ 16 pc ~ 16 Ax and is decreasing 
rapidly.^ Only a very small fraction of the ejecta has made it into 
the dense gas located at (x, y) = ( 0 , 0 ) in the metallicity slice in 
Figure 2. Dilution with metal-free gas has brought the metallicity 
in the densest gas down to below 10 

3.2 Ejecta dispersal trends 

Figure 3 shows the evolution of median radii of Lagrangian tracer 
particles as functions of the original radial mass coordinate in the 
freely-expanding SNR at insertion. Apart from an early retrogres¬ 
sion seen in the innermost shells ~ 20 — 50 kyr after SNR insertion 
that reflects reverse shock passage, all median Lagrangian radii are 
in expansion until ~ 3 Myr. Then the median radii of the envelope 
at mass coordinates >0.3 Mgj begin to stall. Interestingly, the me¬ 
dian radii of the innermost mass coordinates corresponding to the 
core of the explosion continue to expand until the end of the sim¬ 
ulation. This inversion in the median Lagrangian radii in the inner 
10% of the ejecta already begins at ~ 1 Myr. The inversion is a hy¬ 
drodynamic consequence of post-reverse-shock entropy variation 
(Ritter et al. 2015). This is most easily seen in the Sedov-Taylor 
idealization in which pressure tends to a constant value toward the 
center of the remnant whereas density tends to zero steeply with the 
ratio of the radius to the blastwave radius p oc = 

(r/i^sx)^^^ (see, e.g., Landau & Lifshitz 1987, §106; the power 
is 9/2 for non-relativistic electrons with adiabatic index 7 5/3 

and even steeper for relativistic electrons and 7 4/3). Together, 

these scalings imply that specific entropy rises toward the center as 
s oc ln(P/p^) + const oc — ^ ln(r/PsT) + const. 

We can treat the Sedov-Taylor radial coordinate r/PsT oc 
as a function of the fraction m of the swept up mass. 
Then the specific entropy varies with the mass coordinate as 
s oc — Inm + const. Hydrodynamic instability perturbs the 
remnant in the radial direction. If the perturbation is adiabatic 
and the swept up gas remains in approximate pressure equilib¬ 
rium (P = const), the cooling time is proportional to tcooi oc 
P/(Ap^) oc P“^/^e®®/^/A oc (P“^/^/A) where A is the 

usual temperature-dependent cooling coefficient in ergs“^cm^. 
Since A is not a strong function of temperature at SNR interior tem¬ 
peratures > 10^ K, cooling time increases toward small mass coor¬ 
dinates. The gas enriched with the supernova progenitor’s envelope 
cools faster than that enriched with the core. The more rapidly cool¬ 
ing outer mass coordinates form a shell; fragments of this shell lose 
momentum as they run into dense surrounding clouds (here, cosmic 

^ The Jeans length is artificially large thanks to our neglect of H 2 cooling. 
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Figure 4. Left panel: The probability that a particle with initial mass coordinate in the Pop III ejecta becomes incorporated in the central gas clump 

with density ^ 10“^^ g cm“^ (red line) and ^ 10“^^ g cm“^ (green line) that is assumed to be capable of forming Pop II stars at the end of the simulation. 
Right panel: Mollweide projection of the probability of incorporation into the central gas clump with the density threshold of 10“^^ g cm“^ as a function of 
the initial direction of ejection. The probability increases logarithmically from to 10“^ as the color ranges from black to yellow. The radially-averaged 
gas density in the annulus lOpc ^ r ^ 200 pc at SNR insertion is shown in white contours; the contours are logarithmically spaced in density between 
10 “ 24.8 g cjq -3 (dotted line) and 10 “ 23.8 g 0 ^ 1-3 (solid line). 


web streams that have not been photoionized). Shocked gas ele¬ 
ments that do not cool become upward buoyant when the remnant 
comes into pressure equilibrium with the ambient medium. This, 
we believe, is the origin of the inversion in which the innermost 
mass coordinates of the ejecta travel farther than outer ones. 

We turn to enrichment of the densest p ^ 10~^^g cm“^ cen¬ 
tral gas clump that is poised to form second-generation stars. Fig¬ 
ure 4, left panel, shows the return probability, namely, the proba¬ 
bility that a Lagrangian particle becomes incorporated in the dense 
clump, as a function of the radial mass coordinate in the ejecta. The 
probability is markedly reduced for the inner 20% of the ejecta. 
This again is consistent with the hypothesis that the innermost 
ejecta are raised to high entropy, cannot cool, and because of up¬ 
ward buoyancy prefer to reside at larger radii. 

Figure 4, right panel, shows the return probability as a function 
of the initial direction in which the particle was traveling in the 
freely-expanding ejecta. The return probability is very small in the 
darkest colored regions and higher by five orders of magnitude in 
the regions colored bright yellow (the color scale is logarithmic). 
Clearly the variation of the return probability in the ejecta direction 
is much more dramatic than in the ejecta radial mass coordinate. 
The white contours represent the mean density of gas in the annulus 
10 pc ^ r ^ 200 pc as a function of direction. This is the density 
that the ejecta “sees” as it expands into a highly anisotropic gaseous 
environment. The density maxima correspond to filaments of the 
cosmic web and minima correspond to voids. The return probability 
is correlated with the density in the sense that the Lagrangian tracer 
particles shot into directions with higher obstructing densities are 
more likely to become incorporated in the central gas clump. 


3.3 Chemical variation in diffuse and dense gas 

Having identified the gross chemical dispersal trends, we turn to 
the question, to what extent is the long-term gravitational hydrody¬ 
namics of a Pop III SNR able to erase an original chemical stratifi¬ 
cation in the supernova explosion? In Ritter et al. (2015), we found 
that the SNR remained chemically heterogeneous until a fraction 
of the enriched volume was gravitationally-compressed to densi¬ 
ties > 10cm“^. In this dense, gravitationally-collapsing clump, 
gravitational infall excited turbulence in which vorticity was high 
enough to ensure turbulent mixing and local chemical homogeniza¬ 
tion. 

Here, we carry the analysis further by deriving the full statis¬ 
tics of chemical abundance ratios. We divide the ejecta mass co¬ 
ordinate at insertion into 7 equal-mass bins. At any point r we 
compute adaptive Epanechnikov kernel estimates of the local den¬ 
sities pi(r) and Pj(r) of the material originating from bins i and 
j ^ z by choosing the kernel radius Rij (r) such that the combined 
kernel-weighted particle mass within the kernel footprint equals a 
reference mass M, 




peyi^yj ,\rp\<Rij(r) 


1 - 


Ruiry ) 


Trip — M. 


(7) 


Here, yi is the set of Lagrangian tracer particles originating from 
bin i and Vp and rup are, respectively, the position and mass of par¬ 
ticle p. For the reference mass we choose M = 2 x 10“^ M© = 
50 Mpart; this means that on average, 125 particles are located in 
the kernel footprint. Given this adaptive estimate of the kernel ra¬ 
dius Rij{r), we sample space with tracer particles in both bins, 
i and j, and compute kernel-weighted estimates of the densities 
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Figure 5. Distributions of ejecta abundance fractions Xij = Pi/Pj the end of the simulation, where z, j = 1, •••, 7 label the radial mass bins at the time 
of SNR insertion. The histograms were computed with all Lagrangian tracer particles (green lines) and only particles within 50 pc from the point of collapse 
to high densities (red lines). The local ejecta densities were smoothed employing an Epanechnikov kernel with kernel radius adjusted such that the combined, 
kernel-weighted metal mass, as defined in the text, was 2 x 10“^ Mq. Ejecta abundance ratios with | log Xij I > 2 were compressed into the flanking bins. 


Pi{r) and pj (r). We divide the two densities to obtain the ratios 


In Figure 5 we plot the probability density function (PDF) of 
fijO-OgXij) foi" ^^ch pair of ejecta mass coordinate bins. In ad¬ 
dition, we compute the PDF fbo.ijO-OgXij) by restricting kernel 
estimate locations r to Lagrangian tracers that end up within 50 pc 
from the halo center at the end of the simulation. 

The spread in fij decreases with increasing i, showing that 
the outer, higher-initial-velocity mass shells of the ejecta are rel¬ 


atively well-stirred, with the small spread consistent with Poisson 
sampling noise, but the innermost mass shells, however, are poorly 
stirred with each other and the outer shells. The distribution ex¬ 
hibits a heavy tail deficient in inner shell ejecta (corresponding to 
low values of Xij for z = 1 or 2). The halo center distribution 
f 50 ,ij is more narrowly peaked at low i (inner shells) and exhibits a 
similar trend or decreasing spread with increasing ejecta mass co¬ 
ordinate. The displacement of the PDF peak toward low xij shows 
that, consistent with the inner-ejecta-deficiency trend identified in 
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Figure 6. Illustration of possible chemical space effects in Pop II stars resulting from chemical stratification or anisotropy in the Pop III supernova explosion. 
The abundance ratios [Fe/H], [C/Fe], and [Mg/Fe] are compared for three types of artificial model explosions: completely mixed explosions (blue circles), 
spherically symmetric, radially stratified explosions (red triangles), and anisotropic explosions (green squares). The yields of model explosions are calibrated, 
from top to bottom row, to 20, 30, and 40 Mq models of Heger & Woosley (2010) with explosion energies 10^^ erg, mixing parameter 2.51 x 10“^, and 
explosion mechanism S4 (ultra-energetic explosions were selected because spherically symmetric 10^^ erg explosions of Heger & Woosley did not produce 
iron—iron ejection requires either broken spherical symmetry or very high temperatures). For anisotropic explosions we emulated Rayleigh-Taylor fingers 
seen in the Wongwathanarat et al. (2015) core-collapse simulation: each element was shot into a discrete, random set of directions, 4 different directions for 
Fe and 12 directions for Mg. Supernova models were mapped onto tracer particles at insertion and the particles in gas with density ^ 10“^^ g cm“^ at the 
end of the simulation were assumed to be the material from which a Pop II star cluster would form. The purple arrow shows the chemical space displacment 
from fully mixed ejecta to spherically symmetric, radially stratified ejecta. The black + symbols show observed abundance ratios in UFDs (Frebel et al. 2010; 
Ishigaki et al. 2014) and small brown x symbols show those in the Galactic stellar halo (Roederer et al. 2014). 
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Section 3.2, the central 50 pc is indeed deficient in the ejecta from 
the innermost mass coordinates of the explosion. 

3.4 Abundance biases and anomalies in metal-poor stars 

The question that remains is, how should the hydrodynamic trends 
and variation that we have identified in Sections 3.2 and 3.3 mani¬ 
fest in stellar chemical abundance spaces? What is the pure effect of 
explosion geometry—when the yields are themselves held fixed— 
on the enrichment of second-generation stars? Obviously, the single 
explosion that we have simulated cannot answer this general ques¬ 
tion. We can still use the simulation as an indicator of what is the 
plausible magnitude of such effects. For this purpose, we construct 
artificial model explosions: a fully mixed explosion, a spherically- 
symmetric, radially-stratified explosion, and an anisotropic explo¬ 
sion. In search of a suite of template Pop III star explosion models, 
we sought such models in the database of one-dimensional piston- 
triggered explosions of Heger & Woosley (2010). We needed mod¬ 
els that produced at least some iron, but their canonical-energy 
(10^^ erg) explosions did not as all iron group elements were lost 
in the compact remnant. This is of course an artifact of spherical 
symmetry. 

To understand how spherical symmetry breaks down, consider 
the iron core collapse. The fiow onto the proto-neutron star re¬ 
bounds and breaks into an accretion shock. Energy loss to nuclear 
dissociation and neutrino emission stalls the shock, but neutrinos 
from the neutron star heat the shocked, infalling matter and pro¬ 
duce a strong negative radial entropy gradient. The gradient im¬ 
plies vigorous turbulent convection. The largest convective plumes 
reach the standing, already aspherical accretion shock, and in some 
directions push the shock outward, into the envelope of the star. The 
matter in the plumes is dissociated into a-particles and they recom¬ 
bine into ^®Ni as the plumes rise and cool. Nickel is also synthe¬ 
sized in explosive burning of pre-supernova silicon (see, e.g., Ki- 
fonidis et al. 2003). Buoyancy accelerates the under-dense, nickel- 
rich plumes to high velocities. When the plumes cross the inter¬ 
face between the C-0 and He shells, a hydrodynamic process akin 
to the Rayleigh-Taylor instability stretches them into long, narrow 
fingers. Some three-dimensional neutrino-driven, core-collapse ex¬ 
plosions of Wongwathanarat et al. (2013, 2015) and Utrobin et 
al. (2014) that track the explosions until shock breakout develop 
prominent fingers with amplitudes increasing with decreasing mass 
coordinate. The most prominent ^®Ni fingers can outrun the weaker 
i6o, ^^Mg fingers, and both can transport the heavier el¬ 

ements ahead of the isotropically-ejected shell of This ensures 
iron group escape. 

It should be emphasized that the production of nickel fingers 
is sensitive to the precise shock kinematics inside the star; not all 
stellar models in Wongwathanarat et al. (2015) produce fingers. 
When the forward shock crosses the He-H interface (if the pro¬ 
genitor star has retained its hydrogen envelope), it slows down dra¬ 
matically, and this excites a reverse shock propagating backward, 
toward smaller mass coordinates. The reverse shock in turn deceler¬ 
ates any outward-traveling material it crosses. If a finger is slower 
than and is trailing the forward shock, as in the 20 M© blue su¬ 
pergiant model in Wongwathanarat et al., then it collides with the 
reverse shock. This compresses and potentially even disrupts the 
finger. On the other hand, if the finger is faster than and is leading 
the forward shock, as in the 15 M© blue supergiant model, nickel 
is already ahead when the reverse shock forms and is not affected. 
The nickel travels ballistically toward the stellar surface and be¬ 
yond. The finger evolution becomes more complicated in less clear- 


cut cases, such as in both of Wongwathanarat et al.’s red supergiant 
models. 

In view of this multidimensional complexity and the limita¬ 
tions it implies for the robustness of the Heger & Woosley models, 
for our analysis we artificially adopt their ultra-energetic (10^^ erg) 
explosion yield patterns that did produce iron, albeit in small quan¬ 
tities. Specifically, we select three models with progenitor masses 
of 20, 30, and 40 M© and radial stratification as in Figure 7b of 
Heger & Woosley. For anisotropic explosions, we mimic, in a sta¬ 
tistical fashion, the fingering as seen in the Wongwathanarat et al. 
(2015) simulations. Emulating their Figure 3, we take there to be 
4 Fe fingers and 12 Mg fingers. The radial extent of the carbon 
shell was 0.1 — 0.15 Mgj in the ejecta mass coordinate. The an¬ 
gular diameters of the Fe and Mg fingers were 20° and 14°, re¬ 
spectively, and their mass coordinate extents were 0 — 0.3 Mej and 
0.05-0.2 Mej. 

We use the template models to assign chemical identity to La- 
grangian tracer particles based on their location in the SNR at inser¬ 
tion. Then we study chemical abundance ratios [Fe/H], [C/Fe], and 
[Mg/Fe] in the tracer particles that have returned to the halo center 
to become incorporated in the dense ^ 10“^^ gcm“^ clump. In 
anisotropic explosions, we generate several realizations by shoot¬ 
ing fingers in random directions. The resulting abundance ratios are 
shown in Figure 6. All models predict values of [Fe/H] that are too 
low compared to observations, but this is again an artifact of the 
loss of iron to the compact remnant in the spherically-symmetric 
computations of Heger & Woosley. The chemical space displace¬ 
ment from completely mixed to stratified models is 

A[Fe/H] -0.5, 

A[C/Fe] +0.3, 

A [Mg/Fe] +0.1, (9) 

and the displacements in anisotropic models lie in the interval 


(-1.2, -0.8) 

< 

A [Fe/H] 

< (+0.2, +0.3), 


(-0.4, -0.2) 

< 

A [C/Fe] 

< (+0.7,+1.2), 


(-0.7, -0.5) 

< 

r»j 

A [Mg/Fe] 

< (+0.7,+1.1), 

(10) 


where the values in parentheses reflect model-to-model variation. 
The variation in anisotropic models comes close to reproducing the 
scatter in the observed UFDs and the Galactic halo. 


4 DISCUSSION AND CONCLUSIONS 

Previous studies of chemical abundance pattern formation in the 
first Pop II stars (e.g., Karlsson et al. 2013; Ishigaki et al. 2014; 
Cooke & Madau 2014; Kobayashi et al. 2014; Tominaga et al. 
2014) have assumed that individual Pop III supernovae contribute 
their nucleosynthetic yields “monolithically”: the abundances ob¬ 
served in the most metal-poor stars are linear combinations of the 
integral yields of the contributing supernovae. This approach rests 
on the assumption that the supernova ejecta is completely mixed 
either in the explosion, or in the SNR, or perhaps in the larger 
confining gravitational potential well of the host dark matter halo. 
Here and in Ritter et al. (2015), we have provided counterexam¬ 
ples to this assumption, all deriving from the inevitable observation 
that the spatial distribution and kinematic state of the heavier el¬ 
ements returned by the explosion are highly inhomogeneous and 
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anisotropic and that the long-term SNR hydrodynamics can par¬ 
tially preserve such chemical heterogeneity.^ 

One mechanism that we have identified is that the inner ejecta 
cool less efficiently than the outer ejecta due to the inverse relation¬ 
ship between the radius and the amount of specific entropy creation 
in the passage of the reverse shock. This mechanism acts to pre¬ 
serve some of the initial radial stratification of the ejecta in the ma¬ 
terial returning to the halo center to form second-generation stars. 
The other mechanism relates to angular anisotropy of radial decel¬ 
eration as ejecta collide with an anisotropic distribution of dense 
clouds surrounding the explosion. This mechanism preserves an¬ 
gular variation in chemical composition of the ejecta, the variation 
that could arise from, e.g., Rayleigh-Taylor fingering preceding su¬ 
pernova shock breakout. Both mechanisms act to bias stellar en¬ 
richment against the gross yields of the contributing supernovae. 
They can also in themselves contribute to the scatter in abundance 
ratios, potentially explaining a fraction of the observed scatter in 
the most metal-poor stars (most of the scatter is a consequence of 
the stochastic nature of the formation of and nucleosynthesis in Pop 
III and early Pop II stars). Thus, in the initial stage of chemical en¬ 
richment, we predict an irreducible degree of abundance scatter, 
either from anisotropy in a single explosion (this study), or from 
entropy effects in a small multiple of contributing supemovae (Rit¬ 
ter et al. 2015). 

In conclusion, inspired by recent three-dimensional simula¬ 
tions showing extreme inhomogeneity and anisotropy in neutrino- 
powered core-collapse supemovae (Wongwathanarat et al. 2013, 
2015; Utrobin et al. 2014), as well as by the growing evidence 
that the most metal-poor stars exhibit abundance anomalies (Feltz- 
ing et al. 2009; Cohen et al. 2013; Yong et al. 2013) not expected 
for full mixing of supernova yields, we carried out a simulation of 
the long-term transport of nucleosynthetic products from a Pop III 
supernova explosion until ejecta-enriched gas has recollapsed into 
into a dense. Pop II-star-forming clump. Tracing elements back to 
their specific locations in the explosion and mapping onto inho¬ 
mogeneous explosion models, we modeled the scatter in the abun¬ 
dance ratios of the most metal-poor stars. A proof of principle 
rather than a predictive model, our computation shows that SNR 
hydrodynamic effects can contribute to the large scatter in [C/Fe] 
and [Mg/Fe] in the most metal-poor stars. In view of the acceler¬ 
ating stellar chemical abundance data gathering with large spectro¬ 
scopic surveys such as HERMES-GALAH (Zucker et al. 2012; De 
Silva et al. 2015) and Gaia-ESO (Gilmore et al. 2012) as well as 
the recent progress in ab initio supernova modeling (for a review, 
see Mezzacappa et al. 2015, and references therein), Lagrangian 
tagging of the small scale structure of supernova ejecta in simula¬ 
tions of cosmic chemical enrichment and star formation seems like 
a particularly promising tool for future inquiry. 
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